'''
import numpy
#from scipy.optimize import linprog
from scipy.interpolate import CubicSpline
from cpython.ref cimport PyObject
cimport numpy
import sys
cimport cython

cdef public int callPyCubic(a, b, list_c0, list_c1, list_c2, list_c3, int N):
    x = numpy.array(a)
    y = numpy.array(b)
    #print(x)
    #print(y)
    cs = CubicSpline(x, y, bc_type='not-a-knot')
    arr0 = numpy.array(cs.c.tolist()[0])
    arr1 = numpy.array(cs.c.tolist()[1])
    arr2 = numpy.array(cs.c.tolist()[2])
    arr3 = numpy.array(cs.c.tolist()[3])
    for i in range(N-1):
        list_c0[i] = arr0[i]
        list_c1[i] = arr1[i]
        list_c2[i] = arr2[i]
        list_c3[i] = arr3[i]
    return 0
'''

#from typing import Union, Sequence, Tuple, TypeVar, Generic, Optional, List
#from numbers import Number
import numpy
#import functools
import scipy.sparse as sp
import scipy.sparse.linalg as la
from scipy.interpolate import PPoly
from scipy.interpolate import CubicSpline
from cpython.ref cimport PyObject
cimport numpy
import sys
cimport cython
#import array
#from cpython cimport array



def _make_spline(x, y, w, smooth, pcount):
        #pcount > 3
        dx = numpy.diff(x)
        dy = numpy.diff(y)
        dy_dx = dy / dx

        # Create diagonal sparse matrices
        diags_r = numpy.vstack((dx[1:], 2 * (dx[1:] + dx[:-1]), dx[:-1]))
        r = sp.spdiags(diags_r, [-1, 0, 1], pcount - 2, pcount - 2)

        dx_recip = 1. / dx
        diags_qtw = numpy.vstack((dx_recip[:-1], -(dx_recip[1:] + dx_recip[:-1]), dx_recip[1:]))
        diags_sqrw_recip = 1. / numpy.sqrt(w)
        '''
        qtw = (sp.diags(diags_qtw, [0, 1, 2], (pcount - 2, pcount)) @
               sp.diags(diags_sqrw_recip, 0, (pcount, pcount)))
        qtw = qtw @ qtw.T
        '''
        qtw1 = sp.diags(diags_qtw, [0, 1, 2], (pcount - 2, pcount))
        qtw2 = sp.diags(diags_sqrw_recip, 0, (pcount, pcount))
        qtw = numpy.dot(qtw1, qtw2)
        qtw = numpy.dot(qtw, qtw.T)
        p = smooth
        pp = (6. * (1. - p))
        # Solve linear system for the 2nd derivatives
        a = pp * qtw + p * r
        b = numpy.diff(dy_dx).T

        u = la.spsolve(a, b)
        d1 = numpy.diff(numpy.pad(u,(1,1),'constant',constant_values=(0,0))) / dx
        d2 = numpy.diff(numpy.pad(d1,(1,1),'constant',constant_values=(0,0)))

        diags_w_recip = 1. / w
        w = sp.diags(diags_w_recip, 0, (pcount, pcount))
        d3 = pp * w
        d3 = d3 * d2.T
        
        yi = y.T - d3
        pu = p * u
        #pu = numpy.pad(pu,(1,1),'constant',constant_values=(0,0))#(0,0) is equal to 'natural'
        # this is equal to 'not-a-knot'
        acc_start = ((dx[0]+dx[1])*6*pu[0] - dx[0]*6*pu[1]) / dx[1] / 6
        acc_end = ((dx[pcount-2]+dx[pcount-3])*6*pu[pcount-3] - dx[pcount-2]*6*pu[pcount-4])/dx[pcount-3] / 6;
        pu = numpy.pad(pu,(1,1),'constant',constant_values=(acc_start, acc_end))

        c1 = numpy.diff(pu) / dx
        c2 = 3. * pu[:-1]
        c3 = numpy.diff(yi) / dx - dx * (2. * pu[:-1] + pu[1:])
        c4 = yi[:-1]

        #c_shape = (4, pcount - 1)
        #c = numpy.vstack((c1, c2, c3, c4)).reshape(c_shape)
        return c1, c2, c3, c4


cdef public int callPyCsapsF(a, b, w, double p, list_c0, list_c1, list_c2, list_c3, int N) :
    x = numpy.array(a)
    y = numpy.array(b)
    weight = numpy.array(w)
    smooth = p
    c1, c2, c3, c4 = _make_spline(x,y,weight,smooth,N)
    for i in range(N-1):
        list_c0[i] = c1[i]
        list_c1[i] = c2[i]
        list_c2[i] = c3[i]
        list_c3[i] = c4[i]
    return 0


cdef public int callPyCubic(a, b, list_c0, list_c1, list_c2, list_c3, int N):
    x = numpy.array(a)
    y = numpy.array(b)
    #print(x)
    #print(y)
    cs = CubicSpline(x, y, bc_type='not-a-knot')
    arr0 = numpy.array(cs.c.tolist()[0])
    arr1 = numpy.array(cs.c.tolist()[1])
    arr2 = numpy.array(cs.c.tolist()[2])
    arr3 = numpy.array(cs.c.tolist()[3])
    for i in range(N-1):
        list_c0[i] = arr0[i]
        list_c1[i] = arr1[i]
        list_c2[i] = arr2[i]
        list_c3[i] = arr3[i]
    return 0


